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Nonlocal van der Waals density functional: The simpler the better 

Oleg A. Vydrov and Troy Van Voorhis 

Department of Chemistry, Massachusetts Institute of Technology, 
Cambridge, MA, 02139, USA 

(Dated: 27 December 2010) 

We devise a nonlocal correlation energy functional that describes the entire range of dispersion interactions in 
a seamless fashion using only the electron density as input. The new functional is considerably simpler than 
its predecessors of a similar type. The functional has a tractable and robust analytic form that lends itself to 
efficient self-consistent implementation. When paired with an appropriate exchange functional, our nonlocal 
correlation model yields accurate interaction energies of weakly-bound complexes, not only near the energy 
minima but also far from equilibrium. Our model exhibits an outstanding precision at predicting equilibrium 
intcrmonomer separations in van der Waals complexes. It also gives accurate covalent bond lengths and 
atomization energies. Hence the functional proposed in this work is a computationally inexpensive electronic 
structure tool of broad applicability. 



I. INTRODUCTION 

In Kohn-Sham density functional theory (DFT)fi^ 
only one contribution to the ground state energy is not 
known exactly — the correlation energy. Common ap- 
proximations to the correlation energy have the form of 
local or semilocal density functionals^ One of the nu- 
merous limitations of (semi)local functionals is their in- 
herent inability to describe such long-range correlation 
effects as dispersion (van der Waals) interactions. 4 Em- 
pirical dispersion corrections are quite popular and rea- 
sonably successful, but they typically entail a departure 
from pure DFT into the realm of classical force fields*- 
The proper physics of long-range van der Waals inter- 
actions can be captured with fully nonlocal correlation 
functionals. 4 However, the rigor usually comes at the 
cost of an explicit and cumbersome dependence on Kohn- 
Sham orbitals, both occupied and virtual, as exemplified 
by the random-phase approximation and related methods 
(see e.g. Refs.[^and0and references therein). An elegant 
compromise between rigor and computational tractabil- 
ity is achieved in the recently introduced class of nonlocal 
correlation functionals that treat the entire range of dis- 
persion interactions in a general and seamless fashion, 
yet include no explicit orbital dependence and use only 
the electron density as input In this article we devise 
a nonlocal van der Waals functional that is considerably 
simpler yet more accurate than any of its precursors 
We give no derivation for the new functional form, but we 
show that it recovers all the relevant limits and satisfies 
many exact constraints. We suggest a suitable combi- 
nation of exchange and correlation terms that performs 
remarkably well in all our benchmark tests. 



II. FORMALISM 

We write the nonlocal correlation energy as 

Ef = ~ // clrdr' n(r) 1>(r, r') n(r'), (1) 



where n is the total electron density. Building upon the 
insights gained in Refs. Hcl - ll3l we write the correlation 
kernel as 



3e 4 



with 



and similarly 



2m 2 gg'(g + g') 
g = u (v)R 2 + k{t) 



g' = w (r')i? 2 + K (r'). 
In the above equations, R = |r — r'| and 



w (r) = \U) 2 Jt) + 



(2) 



(3) 



(4) 



(5) 



with the local plasma frequency defined via uj 2 
Anne 2 / m and the local band gap given by 



= c- 



Vn(r) 



n(r) 



(6) 



where C is a parameter adjusted to give accurate asymp- 
totic van der Waals Cg coefficients, as elaborated in Sec- 
tion us 

In the R — > oo limit, 



3e 4 



2m 2 Uo (r) Wo (r')[wo(r)+wo(r')]i? 6 ' 



(7) 



so that E~ of Eq. ([T]) has precisely the same lon g-ra nge 
behavior as the nonlocal energy in VV09 of Ref. hll A 
detailed analysis and physical justification of this asymp- 
totic form was given in Ref. [lj. 

For R — > 0, the correlation kernel of Eq. ([2]) behaves 

as 



-A + BR 2 



(8) 
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in accordance with the result of Koide 1 ^ for the R 
asymptotic behavior of the dispersion interaction energy. 
In Eqs. ((3|) and (j4}, we introduced the quantity 



k{t) 



w p (r) 



3b 



fc2(r)' 



(9) 



where Vp — (S^nY^h/m is the local Fermi velocity, 
k s = ^/Sujp/vF is the Thomas-Fermi screening wave vec- 
tor, and b is an adjustable parameter that controls the 
short-range damping of the R~ e asymptote. 
In the uniform density limit, Eq. ([2]) reduces to 



3e 4 
4m 2 



-^R 2 



V 2 



(10) 



and Eq. ([TJ gives the following energy density per elec- 
tron: 



i? 2 $ uni dR 
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32 m 2 v F 



3/4 



32 ao 
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-p, 



(11) 



where ao = Ti 2 /me 2 is the Bohr radius and j3 is a density- 
independent constant. 

It is instructive to rewrite Eq. (1101) in a different form: 



<I>' 



9^3 e 4 
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(12) 



The above equation shows that the R~ 6 asymptote is 
damped at short range on the length scale given by k s R. 
Finally, we define our van der Waals density functional 

as 



E, 



v\ in - E nl + j3N 
dr n(r) 



Jdrn{v) p + ^ /dr'n(r')$(r,r') 



(13) 



where N — J drn(r) is the number of electrons and (3 
is determined by Eq. (1TT1) . By construction, E^ V1 ° of 
Eq. ()13|) vanishes in the uniform density limit. This is 
a useful property: we can pair E^ V1 ° with an existing 
exchange-correlation (XC) functional without affecting 
the description of the uniform electron gas. 

It is a nontrivial task to determine how Eq. (fl"3| be- 
haves in the slowly varying density limit. At any rate, 
it has been argued that recovery of the density gradient 
expansion for correlation energy is of little relevance for 
real systems.— 

The van der Waals functional of Eq. (fl~3|) is very simple 
and easily implementable. The double integral in Eq. (fTJ) 
is in practice evaluated as a double sum over a numer- 
ical grid. For R — > 0, the kernel of Eq. ([2]) goes to a 
constant, avoiding any complications with the numerical 
integration. In the inner integration loop, only simple 



arithmetic operations need to be performed. Hence this 
new functional is computationally less expensive than the 
nonlocal functionals of Refs. la-fllT 

The previous version of the nonlocal van der Waals 
functional, VV09 proposed by us in Ref. [llj, was more 
computationally demanding because it used a rather 
elaborate damping function. Evaluation of the nonlocal 
correlation energy in VV09 required computing a square 
root, an exponent, and an error function for each pair of 
grid points in the double sum. The energy derivatives, 
needed for the self-consistent implementation of VV09, 
were also quite complicated. 16 With the new VV10 model 
of Eq. (Tl3l) . these energy derivatives are much simpli- 
fied. The self-consistent implementation of VV10 within 
a Gaussian basis set code is described in detail in the 
next section. 



III. IMPLEMENTATION 

For the sake of brevity, in this section we switch to 
atomic units (h — m — e = 1). In these units, n of 
Eq. ([9J is expressed rather simply as 



k(t) 



3tt 



n(r) 



-l 1/6 



9tt 



(14) 



Our implementation of VV10 is very similar to the 
implementation of its predecessor VV09, described in 
Ref. O We express the electron density via a Gaussian 
basis set: 



(15) 



where Xfj. an d Xv are basis functions, and P^ v are the den- 
sity matrix elements. For the self-consistent treatment, 
we need to find the derivatives of E^ V1 ° with respect to 



d£ c vvl ° 
dP, 



SE VV1Q 

on(r) 



(16) 



To that end, we employ the standard formalism 1 ^ devel- 
oped for semilocal XC functionals: 



d£ c vvl ° 
dP„„ 



dr 



F„x»Xu+2F-yVn-V(xnX») ■ (17) 



We denote 7 = |Vn| 2 for convenience. F n and F 7 in 
Eq. p7|) arc computed as 

F n (r) = j3 + J dr'n(r')$ 

+ n(r) 

and 



, s TT , \ Ouj , . TT ,, , 

a^ (r) (r) + TfrT (r) (r) 



(18) 



F 7 (r)=n(r)^(r)W(r), 



(19) 
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where 



U(r) = -J dr'n(r')$ 



and 



W(r) = - J dr' ra(r') |r - r'f $ 
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9 + .9' 



(20) 



(21) 



In the above equations, $, g, and g' are assumed to be 
functions of both r and r' . 

Within an atom-centered basis set implementation, the 
gradient of E^ yw with respect to the displacement of 
nucleus A has three contributions: 



A A A 

SGBF + Swcights + Sgrid ■ 



(22) 



Sgbf denotes the contribution of the Gaussian basis 
functions. This term can be evaluated by plugging F n 
and F 1 into Eq. (9) of Ref. [13 instead of df/dn and 

df/d 1 . 

The last two terms in Eq. (|2"2"|) are due to the use 
of atom-centered numerical integration quadratures. We 
employ the atomic partitioning scheme of Becke/^ which 
separates the molecular integral into atomic contribu- 
tions: 



£c vvl0 = EE^"(^) 



A i£A 



\ EE WB i n ( r Bj) Hn- Ai ,r Bj ) 



b jeB 



(23) 



where wai and WBj are the quadrature weights, and the 
grid points va% are given by r a% — Ra + where Ra 
is the position of nucleus A, with the defining a one- 
center integration grid. The quadrature weights depend 
on the nuclear configuration and hence have a nonzero 
gradient with respect to nuclear displacements: 



Swcights 



EE (^AW B i)n{r B i) 



B i£B 



P 



(24) 



+ E E w °3 n ( r Cj) ®( r Bi, r Cj ) 
c jec J 

The weight derivatives V awbi can be found in Ref. 

The last term in Eq. (|22|) arises because each one- 
center quadrature grid moves together with its parent 
nucleus and the nonlocal correlation kernel $ depends 
explicitly on the distance between the grid points. The 

Sgrid 



v ' term can be computed as: 



i£A 

X E ^ WBj n(rBj)Q(TAh *Bj) {YAi - TBj) , 
B=£Aj£B 

(25) 



where 
Q(r,r') 



2$ 



w (r) wo(r') w (r) + w (r') 



9 9' 9 + 9' 

(26) 

Analytic gradients, computed via Eq. (|2"2"1) . enable us 
to perform structural optimizations routinely and effi- 
ciently. 



IV. ADJUSTMENTS 

^vvio Q f p^|) can in principle be paired with 
nearly any existing XC functional. However, to avoid 
double-counting as much as possible, it is preferable to 
combine E^ YW with a functional that gives no significant 
binding in van der Waals complexes. Since predictions of 
XC functionals for weakly interacting systems can differ 
drastically, the parameter b, controlling the short-range 
behavior of E^ V1Q , has to be adjusted for a particular 
parent XC approximation. 

The Perdew-Wang 86 (PW86) exchange functionate 
has been shown to describe the repulsive parts of van der 
Waals potentials rather well^Sr— and a refitted version 
of PW86 was recently proposed^ We denote this 're- 
fitted PW86' as rPW86 here. In the rest of this article 
we will consider one particularly apt combination of the 
exchange and correlation terms: 



pVVlO prPW86 



£PBE + £ ,VV10 ) (27) 

where E^ BE is the semilocal correlation energy functional 
in the generalized gradient approximation of Perdew, 
Burke, and Ernzcrhof (PBE) M- Hereafter, we will often 
refer to the XC functional of Eq. (gZJ) simply as VV10. 

There are two adjustable parameters in our method, 
but only C in Eq. ([6]) affects the asymptotic dispersion 
Cq coefficients. We fit C to minimize the average error 
for the benchmark set of 54 Cq coefficients compiled in 
Tables II and III of Ref. [IH Using self-consistent electron 
densities from rPW86-PBE (i.e. using E xc = ££ PW86 + 
£ C PBE ), we find the optimal value of C = 0.0093, which 
gives the mean absolute percentage error (MAPE) of 14% 
for the whole test set of 54 species. The largest errors in 
Cq occur for metal atoms. When all the metal atoms (Li, 
Na, Mg, Al, Zn, Ga) are removed from the benchmark 
set and only the remaining 48 species are considered, the 
MAPE in C 6 drops to 9%. Note that a slightly smaller 
value of the parameter C was used in Refs. [U [H, and 
[iH since the source of electron densities was different!^ 

Another adjustable parameter b, introduced in Eq. (0), 
controls the short range damping of Using -E^X 10 of 
Eq. (|2"T|) . we fitted b on the training set of 22 binding 
energies (the S22 set of Ref. l25h . with the computational 
details described in Section [V] and the results given in 
Section IVI Al We obtained the best fit with b = 5.9. 
Correspondingly, in Eq. (|13p we use 



0=^ 
' 32 



5.9 2 



3/4 



= 0.00497 
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in the Hartree atomic units. 

Although the rPW86 exchange functional appears to 
be a very good match for our correlation model, this 
choice is certainly not unique. In the Appendix we de- 
scribe an alternative choice of exchange which is nearly 
equally suitable for this purpose. 



V. COMPUTATIONAL DETAILS 

VV10 has been implemented sclf-consistently into the 
Q-Chem software package.— In our assessment of VV10, 
we compare its performance to another recent van der 
Waals functional of a similar type — vdW-DF2 of Ref. |9|. 
Since vdW-DF2 is just a reparameterization of an ear- 
lier functional, vdW-DF of Ref. we use the same 
implementation^ for both of them. The full XC energy 
in vdW-DF2 is defined as 



pvdW-DF2 _ prPW86 , pLDA 

XC X ' c 



E. 



vdW-DF2 
il ) 



(28) 



where i?^ DA is the correlation energy in the local density 
approximation, for which we use the parameterization of 
Perdew and Wang.— 

All reported calculations were performed self- 
consistently. For weakly bound systems, the interaction 
energies were counterpoise corrected. The unpruned 
Euler-Maclaurin-Lebedev (99,590) quadrature grid was 
used to evaluate all local and semilocal contributions 

(£rPW86 ; E PBE t and ^LDA) and the (75^) grid 

was used for the nonlocal components. Deviations of 
computed binding energies from the reference values 
are analyzed with the help of mean errors (ME), mean 
absolute errors (MAE), and mean absolute percentage 
errors (MAPE). In all tables, binding energies are re- 
ported as positive values. Hence a negative ME indicates 
undcrbinding while a positive ME means overbinding. 

In our recent article^ assessing the performance of 
VV09, we used the same benchmark sets of weakly bound 
systems and the same computational details. Therefore, 
the results reported in Ref. [TH can be directly compared 
to the results obtained with the new functionals, pre- 
sented in Section LYU below. 



TABLE I. Binding energies (in kcal/mol) for the S22 test 
set. All calculations were performed self-consistently with the 
aug-cc-pVTZ basis set and counterpoise corrected. Molecular 
structures are from Ref. [25] and reference binding energies are 
from Ref. [2^ for all systems except Adenine-Thymine com- 
plexes, for which the numbers from Ref. were used. 



Complex (symmetry) 


Ref. 


VV10 


vdW-DF2 


Dispersion-bound complexes (8) 








CH 4 dimer (L> 3 d) 


0.53 


0.50 


0.68 


C2H4 dimer (£>2d) 


1.48 


1.42 


1.32 


Benzene-CBU (C3) 


1.45 


1.45 


1.29 


Benzene dimer (C2h)^ 


2.66 


2.71 


2.15 


Pyrazine dimer (C s ) 
Uracil dimer (C 2 £ 


4.26 


4.02 


3.30 


9.78 


9.70 


8.76 


Indole-Benzene (Ci)_ 


4.52 


4.54 


3.44 


Adenine-Thymine (Ci)^ 


11.66 


11.42 


9.58 


Mixed complexes (7) 








C2H4-C2H2 (C 2 v) 


1.50 


1.68 


1.53 


Benzene-tbO (C s ) 


3.28 


3.31 


2.80 


Benzene-NH 3 (C s ) 


2.32 


2.28 


1.99 


Benzene-HCN (C B ) 


4.54 


4.30 


3.55 


Benzene dimer (C^v)^, 


2.72 


2.54 


2.06 


Indole-Benzene (C s ^ 


5.63 


5.27 


4.20 


Phenol dimer (Ci) 


7.10 


6.99 


5.97 


Hydrogen-bonded complexes (7) 








NH 3 dimer (C 2h ) 


3.15 


3.43 


2.97 


H 2 dimer (C s ) 


5.00 


5.50 


4.78 


Formic acid dimer (C21O 


18.75 


19.96 


16.77 


Formamide dimer (C2h) 
Uracil dimer (C2h)^ 


16.06 


16.71 


14.43 


20.64 


21.10 


18.69 


2-pyridone-2-aminopyridine (Ci) 
Adenine-Thymine WC (d)^ 


16.94 


18.05 


15.37 


16.74 


17.42 


14.74 



a 'Parallel-displaced' configuration. 
b Stacked configuration. 
c T-shaped configuration. 
d Planar configuration. 



VI. RESULTS FOR WEAKLY INTERACTING SYSTEMS 
A. The S22 benchmark set 

Binding energies for the S22 test set, computed with 
vdW-DF2 and VV10 at the geometries of Ref. [H, are 
given in Table [I] The error statistics are summarized in 
Table [TTJ Updated reference values of binding energies 
for the S22 set were recently obtained in Refs. and 
I30L For our benchmark set, we picked the values that we 
deemed the most converged with respect to the basis set 
size: we took the reference values for most of the systems 



TABLE II. Summary of deviations from the reference values 
of the binding energies reported in Table [I] ME and MAE are 
in kcal/mol, MAPE is in percents. 





VV10 


vdW-DF2 


Dispersion-bound complexes (8) 




ME 


-0.07 


-0.73 


MAE 


0.09 


0.76 


MAPE (%) 


2.6 


18.0 


Mixed complexes 


(7) 




ME 


-0.10 


-0.71 


MAE 


0.16 


0.72 


MAPE (%) 


4.8 


16.8 


Hydrogen-bonded complexes (7) 




ME 


0.70 


-1.36 


MAE 


0.70 


1.36 


MAPE (%) 


6.1 


8.8 


Total (22) 






ME 


0.16 


-0.92 


MAE 


0.31 


0.94 


MAPE (%) 


4.4 


14.7 
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from Ref. l29l . except for the Adenine-Thymine complexes 
(in both stacked and planar configurations), for which we 
used the numbers from Ref. [3(J 

As Tables Q] and |TT] show, VV10 performs remarkably 
well for all dispersion-bound and mixed complexes, but 
yields somewhat larger errors for hydrogen-bonded sys- 
tems. The parameter b inside £J^ V1 ° was fitted to the 
binding energies of the S22 set, but it is by no means 
trivial to achieve such accuracy by adjusting just one 
global parameter. 

vdW-DF2 overbinds some of the smaller systems, 
specifically (CH^ and C2H4-C2H2, but underbinds all 
other complexes in the S22 set, when the geometries from 
Ref. [ID are used. This under binding tendency intensifies 
as the monomer sizes increase. A negative ME for the 
binding energies of the S22 set was also found in Ref. 0. 
However, intermonomer separations were optimized in 
Ref. 0, yielding larger binding energies, as compared to 
our results for the fixed^ geometries. 

Binding energies calculated for a standard set of near- 
equilibrium geometries certainly do not tell the full story 
about the performance of a functional. For a more 
comprehensive assessment, it is necessary to determine 
whether equilibrium intermonomer separations are accu- 
rately located and whether reasonable interaction ener- 
gies are predicted far from equilibrium. To this end, we 
analyze binding energy curves for several complexes. 

B. Interaction energy curves 

We computed interaction energy curves for six weakly- 
bound systems using VV10 and vdW-DF2. Figure Q] 
shows the results for the argon and krypton dimers, Fig- 



ure [21 for the methane dimer (D^ symmetry) and the 
benzene-methane complex (C3 V ), and Figure [3] for ben- 
zene dimers in two configurations: T-shaped (C?2v) an d 
stacked sandwich-shaped (D^). We applied counter- 
poise corrections at all intermonomer distances. For Ar2 
and Kr2, the reference values were taken from Ref. [3l|, 
and for the other complexes from Ref. [32|. For molecular 
complexes, we used the same fixed monomer geometries 
as in Ref. |32|. 

For all six complexes, VV10 reproduces the equilib- 
rium intermonomer separations with remarkable preci- 
sion. Interaction energies yielded by VV10 agree quite 
well with the reference values at all separations. The 
well depths are predicted with reasonable accuracy: 
VV10 slightly underbinds (CH^ and the T-shaped 
(CeH 6 )2, but it slightly overbinds Ar 2 and Kr 2 and more 
strongly overbinds the sandwich-shaped (CeHe)2- For 
the benzene-methane complex, the deviations of W10 
from the reference curve are within the uncertainties of 
the reference^ interaction energies at nearly all sepa- 
rations. The long-range falloff of the potential energy 
curves is reproduced very well by VV10 in nearly all 
cases. 

vdW-DF2 tends to predict the energy minima at some- 
what too large separations. It significantly overbinds the 
complexes of small monomers: Ar2, Kr2, and (CIL^. 
For all systems in Figures [1] [2] and|3l except Ar2, the re- 
pulsive walls given by vdW-DF2 are too steep and shifted 
towards larger distances as compared to the reference 
curves. As was shown in Ref. [IH, vdW-DF2 strongly un- 
derestimates asymptotic Cq coefficients. As a result, for 
systems where the asymptotic interactions are dominated 
by dispersion, vdW-DF2 will yield too shallow potential 
energy curves at long range. This underestimation is vis- 
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R(A) R(A) 



FIG. 2. Interaction energy curves for (a) the methane dimer and (b) the benzene-methane complex. R is the distance 
between the centers of mass of the monomers. Reference values are from Ref. I32I VV10 and vdW-DF2 results were obtained 
self-consistently with the aug-cc-pVTZ basis set. 




4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 



R(A) R(A) 

FIG. 3. Interaction energy curves for (a) the T-shaped and (b) the stacked sandwich-shaped benzene dimers. R is the distance 
between the centers of mass of the monomers. Reference values are from Ref. I3H VV10 and vdW-DF2 results were obtained 
self-consistently with the aug-cc-pVTZ basis set. 



ible at large separations for Ar2, Kr2, and the T-shaped 
(CqH$)2- The effects of the poor vdW-DF2 asymptotics 
are expected to be more noticeable at very large inter- 
monomer distances or for very large monomers. Such 
cases are not significantly represented in this study. 

As Figure H(b) shows, both VV10 and vdW-DF2 
significantly overestimate the equilibrium binding en- 
ergy in the sandwich-shaped benzene dimer. This may 



be partially caused by the lack of Axilrod- Teller-type 
three-body interactions in the functionals of the form 
of Eq. (JlJ. Another possible source of errors is the 
rPW86 exchange functional. It was found in Ref. that 
PW86 exchange is more repulsive than Hartree-Fock for 
the T-shaped benzene dimer, but less repulsive for the 
sandwich-shaped configuration. 
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TABLE III. Total energies of atoms in Hartree atomic units. 
The exact values are from Refs. [33| and 0. 





Exact 


PBE 


rPW86-PBE 


VV10 


H 


-0.500 


-0.500 


-0.509 


-0.505 


He 


-2.904 


-2.893 


-2.926 


-2.916 


Li 


-7.478 


-7.462 


-7.517 


-7.503 


Be 


-14.667 


-14.630 


-14.709 


-14.692 


B 


-24.654 


-24.612 


-24.723 


-24.700 


C 


-37.845 


-37.799 


-37.941 


-37.914 


N 


-54.589 


-54.536 


-54.708 


-54.677 





-75.067 


-75.015 


-75.230 


-75.194 


F 


-99.734 


-99.676 


-99.932 


-99.891 


Ne 


-128.938 


-128.866 


-129.159 


-129.114 


Na. 


-162.255 


-162.172 


-162.502 


-162.453 


Mg 


-200.053 


-199.955 


-200.325 


-200.271 


Al 


-242.346 


-242.236 


-242.645 


-242.588 


Si 


-289.359 


-289.234 


-289.682 


-289.621 


P 


-341.259 


-341.115 


-341.602 


-341.536 


S 


-398.110 


-397.952 


-398.481 


-398.411 


CI 


-460.148 


-459.974 


-460.544 


-460.470 


Ar 


-527.540 


-527.346 


-527.955 


-527.876 


ME/? 




0.008 


-0.019 


-0.014 



a Mean error per electron. 



VII. RESULTS FOR ATOMS AND COVALENT BONDS 

We have shown that VV10 describes weakly-bound sys- 
tems quite well. In this section, we show that VV10 is a 
broadly applicable method that also treats covalent and 
ionic bonds well. 



A. Total energies of atoms 

In Table Mil we report the total energies of atoms 
up to Ar, computed with rPW86-PBE (that is rPW86 
exchange^ with PBE correlation^,) as well as VV10, 
which is rPW86-PBE with the addition of E^ V1 °, as 
shown in Eq. ([27| . We also include the popular semilocal 
XC functional PBE<22 We used the aug-pc-4 basis set^ 
with all beyond-/ polarization functions removed. This 
basis set is expected to yield energies close to the basis 
set limit. As the benchmark, we use the accurate non- 
relativistic values of atomic energies from Refs. [33| and 

Among the three considered functionals, PBE gives the 
most accurate atomic energies. PBE total energies are 
systematically above the exact values, except for the H 
atom, for which PBE gives a nearly exact energy. On 
the contrary, rPW86-PBE and VV10 yield total energies 
that are systematically too low. 

We find £VVio of Eq to be positive not only for 
atoms, but for all systems considered in this study. We 
have no proof that E^ V1 ° is always strictly nonnegative, 
i.e. that j3N > —Efr, but we have found no exceptions 
so far. The positive contribution of E^ V1 ° brings the 
rPW86-PBE total energies in better agreement with the 



TABLE IV. Atomization energies for the AE6 set of Ref. H 
computed with the aug-cc-pVTZ basis set. All values are in 
kcal/mol. 





Ref. 


PBE 


rPW86-PRE 


VV10 


s 2 


101.7 


112.5 


106.0 


108.5 


SiO 


192.1 


193.4 


190.1 


191.8 


SiH 4 


322.4 


311.6 


308.2 


310.2 


glyoxal 


633.4 


661.7 


643.9 


650.1 


propyne 


704.8 


720.1 


706.1 


711.5 


cyclobutane 


1149.0 


1166.1 


1137.0 


1148.8 


ME 




10.4 


-2.0 


2.9 


MAE 




14.0 


7.4 


7.2 



exact values, as Table Hill shows. 

We did not include vdW-DF2 in Table [TTT1 because this 
functional is not defined for open-shell systems. For the 
five closed-shell atoms in the set of Table IIII1 vdW-DF2 
yields the ME of -0.048 hartree per electron. vdW-DF2 
total energies are much too low primarily because this 
functional includes £^ DA as a constituent and LDA is 
known to severely overestimate the magnitudes of corre- 
lation energies in atoms. 



B. Atomization energies 

We computed atomization energies for the AE6 set de- 
veloped by Lynch and Truhlar.— This set includes only 
six molecules, but it is quite diverse and was constructed 
to be representative, that is to reproduce mean errors in 
much larger sets. We used the molecular geometries sug- 
gested by the authors of the set^ The results and the 
mean errors are given in Table IIV1 

rPW86-PBE, which is the backbone of VV10, yields 
considerably more accurate atomization energies as com- 
pared to the popular PBE functional. The addition of 
the nonlocal correlation term £^ V1 ° to rPW86-PBE in- 
creases the atomization energies in all cases, as the VV10 
results in Table IIVI show. In other words, the addition 
of E^? vw makes covalent bonds somewhat stronger. As 
a result, VV10 yields smaller errors than rPW86-PBE 
for three molecules in the AE6 set, but larger errors for 
the other three molecules. On average, the performance 
of VV10 for the AE6 set is quite similar to its parent 
functional rPW86-PBE. 



C. Bond lengths 

In Table [V] we have compiled a test set of 20 small 
closed-shell molecules for which experimental 37 equilib- 
rium bond lengths (r e ) are known. The test set con- 
sists mostly of diatomics, but includes several polyatomic 
molecules of high symmetry whose geometry is com- 
pletely determined by a single bond length. Most of the 
molecules in the set are covalently bound, except for the 
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TABLE V. Equi 


librium bond lenj 


jths (r P ) 


computed using the aug 


;-cc-pVTZ basis set. 


All values are in A. Experimental data 


are from Ref. [37 
















Expt. 


PBE 


rPW86-LDA 


vdW-DF2 


rPW86-PBE 


VV10 


H 2 


0.741 


0.751 


0.734 


0.736 


0.745 


0.745 


LiF 


1.564 


1.583 


1.592 


1.592 


1.586 


1.585 


LiCl 


2.021 


2.028 


2.041 


2.041 


2.031 


2.031 


Li 2 


2.673 


2.726 


2.663 


2.670 


2.699 


2.700 


CH 4 


1.087 


1.096 


1.088 


1.090 


1.092 


1.093 


CO 


1.128 


1.137 


1.136 


1.136 


1.137 


1.137 


C0 2 


1.160 


1.172 


1.174 


1.173 


1.172 


1.172 


cs 2 


1.553 


1.566 


1.573 


1.571 


1.568 


1.567 


N 2 


1.098 


1.103 


1.101 


1.101 


1.102 


1.102 


HF 


0.917 


0.932 


0.934 


0.935 


0.934 


0.934 


F 2 


1.412 


1.413 


1.459 


1.453 


1.433 


1.432 


NaCl 


2.361 


2.385 


2.410 


2.401 


2.390 


2.385 


NaBr 


2.502 


2.528 


2.560 


2.547 


2.538 


2.532 


Na 2 


3.079 


3.079 


3.006 


2.987 


3.035 


3.024 


SiO 


1.510 


1.536 


1.543 


1.540 


1.539 


1.538 


P 2 


1.893 


1.911 


1.922 


1.918 


1.914 


1.913 


HC1 


1.275 


1.291 


1.289 


1.290 


1.290 


1.290 


Cl a 


1.988 


2.019 


2.085 


2.073 


2.040 


2.037 


HBr 


1.415 


1.433 


1.436 


1.437 


1.435 


1.435 


Br 2 


2.281 


2.311 


2.396 


2.380 


2.343 


2.336 


ME 




0.017 


0.024 


0.021 


0.018 


0.017 


MAE 




0.017 


0.033 


0.031 


0.023 


0.022 



alkali metal halides where the bonding is ionic. 

Table [V] reports the results of geometry optimizations 
carried out using analytic gradients and the aug-cc-pVTZ 
basis set. We tested vdW-DF2 and VV10 alongside their 
semilocal parent functionals, rPW86-LDA and rPW86- 
PBE respectively. We also included PBE for the sake 
of comparison. For most molecules in Table [V] vdW- 
DF2 yields nearly the same bond lengths as rPW86-LDA, 
and VV10 nearly the same as rPW86-PBE. Hence the 
addition of a nonlocal correlation term has an almost 
negligible effect on covalent and ionic bond lengths. The 
largest change occurs for Na2: the Na-Na bond length 
given by VV10 is 0.011 A shorter than with rPW86-PBE. 
This is not surprising since the Na-Na bond is weaker 
than most covalent bonds. On average for the test set of 
Table [V] the performance of VV10 is on par with PBE, 
while vdW-DF2 is somewhat less accurate. 



VIII. CONCLUSIONS 

The VV10 correlation model proposed in this work be- 
longs to the family^— of nonlocal van der Waals density 
functionals defined in terms of the electron density alone 
and using no orbital input. These functional^— are 
general and seamless: they require neither splitting the 
system into interacting fragments nor any kind of atomic 
partitioning. They treat inter- and intra-molecular dis- 
persion interactions on equal footing. 

VV10 has the same long-range behavior as its precur- 
sor VV09, but the damping mechanism of dispersion in- 
teractions at short range is greatly simplified in VV10. 



This simplification not only makes the functional more 
efficient and computationally tractable, but it also leads 
to improved overall accuracy. 

An essential aspect of the VV10 formalism is the ad- 
ditional flexibility introduced with the help of an ad- 
justable parameter b which controls the short range be- 
havior of the nonlocal correlation energy. When E^ YW 
of Eq. (|T3"|) is added as a correction to an existing XC 
functional, b is adjusted to attain a balanced merging of 
interaction energy contributions at short and intermedi- 
ate ranges. A particularly successful functional is con- 
structed by adding E™° with b = 5.9 to rPW86-PBE, 
as shown in Eq. (|27|) . Another parameter C — 0.0093 
inside E^ YW ensures that accurate asymptotic Cq coef- 
ficients are obtained^ 

As our benchmark tests clearly demonstrate, the func- 
tional of Eq. (|2"T)) is a broadly accurate electronic struc- 
ture tool: its outstanding predictive power for weakly- 
bound systems is complemented by its good description 
of covalent bonds. 
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Appendix: VV10 with long-range corrected exchange 

The XC functional of Eq. (|2"T1) includes the semilo- 
cal rPW86 exchange term as a constituent. Semilo- 
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TABLE VI. Errors of LC-VV10 for several test sets. LC- 
VV10 is defined by Eq. (IXTT) with uj = 0.45, C = 0.0089, and 
6 = 6.3. 





ME 


MAE 


MAPE (%) 


S22 subsets (kcal/mol) 








Dispersion-bound 


0.01 


0.17 


7.1 


Mixed 


0.09 


0.16 


4.0 


Hydrogen-bonded 


0.17 


0.31 


2.3 


Full S22 set 


0.09 


0.21 


4.6 


Other test sets 








Energies of atoms (a.u.) 


0.008^ 


0.008^ 




AE6 (kcal/mol) 


0.5 


6.9 




Bond lengths (A) 


-0.006 


0.014 




a Errors per electron. 



cal exchange functionals are known to suffer from self- 
interaction error (SIE), which causes such artifacts as 
poor description of charge transfer complexes and tran- 
sition states of chemical reactions. Long-range corrected 
(LC) hybrid exchange functionals have proven to be very 
effective at minimizing the SIEi 38 i 39 LC hybrids have also 
been shown to adequately describe the repulsive parts of 
van der Waals potentials^ We test one particular LC hy- 
brid for its suitability as a counterpart for the VV10 cor- 
relation model — the LC-wPBE functional] 38 ' 41 The ex- 
change component in LC-wPBE is defined as a sum of two 
parts: the long-range Hartree-Fock (LR-HF) part and 
the short-range PBE (SR-PBE) part, for which we use 
the parameterization of Ref. |41| . For the range-separation 
parameter ui, we use the value of 0.45 bohr -1 , as sug- 
gested in Ref. Hi] (this variant was termed LC-cjPBE08 
therein). The full XC functional, which we denote LC- 
W10, has the form 



E. 



LC-VV10 



El 



SR-PBE 



.£LR-HF + £PBE +£ ™ (A1) 

The two parameters inside E^ vw are adjusted using the 
same procedure as described in Section IIVI As we pre- 
viously found) 11 ' 13 ' 24 C = 0.0089 gives accurate Cq co- 
efficients at LC-wPBE electron densities. We fit the pa- 
rameter b to the binding energies of the S22 set using the 
functional of Eq. (IA.1[) and obtain the optimal value of 
6 = 6.3. 

In Table IVII we summarize the errors of LC-VV10 for 
the test sets used previously in this article. As com- 
pared to the VV10 model with the rPW86 exchange, LC- 
VV10 performs somewhat better for binding energies of 
hydrogen bonded complexes, for total energies of atoms, 
and for bond lengths. However, no improvement is ob- 
served for atomization energies, while binding energies of 
dispersion-bound complexes are slightly worsened on av- 
erage. Therefore, we suggest using LC-VV10 only when 
SIE is an issue. In most other cases the VV10 model 
of Eq. (|2"T1) is preferable since the semilocal rPW86 ex- 
change is computationally cheaper and much more easy 
to implement than an LC hybrid. 

It is noteworthy that our model can achieve good per- 
formance using a pre-existing and unmodified exchange 



functional (rPW86 or LC-ojPBE) and adjusting only the 
parameters in the nonlocal correlation. It might be pos- 
sible to further improve the performance by tailoring an 
empirical exchange functional for our specific purpose, in 
the vein of the recent works of Refs. 1421444 However, we 
prefer to keep empirical fitting to a minimum. 
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